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We use a nonequilibrium implementation of extended dynamical mean field theory to study the effect of 
dynamical screening in photo-excited Mott insulators. The insertion of doublons and holes adds low-energy 
screening modes and leads to a reduction of the Mott gap. The coupling to low-energy bosonic modes further¬ 
more opens new relaxation channels and significantly speeds up the thermalization process. We also investigate 
the effect of the energy distribution of the doped carriers on the screening. 


I. INTRODUCTION 

Driving a material out of equilibrium by a strong laser 
pulse can provide new insights into correlation phenomena?® 
and even induce transitions into nonthermal “dark” state sP 
A conceptually rather simple example of a pulse induced 
nonequilibrium phase transition is the photo-doping of a Mott 
insulator. 7 ® In these experiments, a pulse with frequency 
larger than the Mott gap produces doublon-hole pairs and 
hence results in a nonthermal conducting state. While the met¬ 
allization happens on a femto-second timescale, the relaxation 
back to a Mott insulating state can take several picoseconds . 7 

On the theoretical side, different aspects of the photo¬ 
doping and thermalization process have been recently investi¬ 
gated. In a system with sufficiently large gap, the relaxation 
proceeds in two stages. In the first stage, the photo-doped car¬ 
riers “thermalize” within the Hubbard band due to electron- 
electron scattering, or loose their kinetic energy due to scat¬ 
tering with phonons 11 12 or spins. 'SSI Scattering processes 
which change the number of doublon-hole pairs, nece ssary 
for thermalization, depend exponentially on the gap size.® I7 i 
In small-gap insulators, the initial kinetic energy of the photo- 
doped carriers may be sufficient for “impact ionization,’® 
which introduces an additional timescale into the problem. 

An important aspect, which has not been considered in 
these previous studies based on the nonequilibrium dynam¬ 
ical mean field approximation , 19 is the effect of screening 
from longer-ranged Coulomb interactions. The interaction pa¬ 
rameters used within an (extended) Hubbard-model type de¬ 
scription are partially screened interactions, which can be ob¬ 
tained from the fully screened interaction by removing the 
screening processes withing the subspace of the low-energy 
model . 20 The low-energy screening is however very different 
in a Mott insulator and in a metal. If mobile carriers are in¬ 
serted into a Mott insulator by photo-doping, the screening of 
the Coulomb interaction will change, and this should have a 
noticeable effect on the nature of the photo-doped state and 
on the relaxation processes. Other effects of the longer-range 
Coulomb interaction on the non-equilibrium dynamics have 
been studied in a one dimensional chain using the time depen¬ 
dent Lanczos method. There, a local enhancement of charge 
(spin) order 21 22 and an appearance of in-gap states 23 can be 
observed. 

A method which captures the screening from long-range 
Coulomb interactions is extended dynamical mean-field the¬ 


ory (EDMFT). While this formalism has been developed more 
than a decade ago,®® an accurate numerical implementa¬ 
tion has only recently become feasible . 27 29 Applications to 
the Hubbard model with on-site and inter-site interactions 
have clarified the phase diagram and the dominant low-energy 
screening modes in the insulating and metallic phase.® Here, 
we extend EDMFT to the nonequilibrium domain by imple¬ 
menting the scheme on a Kadanoff-Baym contour. This al¬ 
lows us to study the dynamical screening effect after a photo¬ 
doping excitation in real time. We show that the screening, 
and the associated possibility to emit and absorb plasmons, 
influences the thermalization process in significant ways and 
that the doping-induced screening of the Coulomb interaction 
has an effect on the gap size of photo-doped Mott insulators. 
We also show evidence for nontrivial transient states induced 
by time-dependent changes of the screening environment. 

This paper is organized as follows: In section II we discuss 
the nonequilibrium generalization of EDMFT and its imple¬ 
mentation using an impurity solver which combines the non¬ 
crossing approximation with a weak-coupling expansion for 
the retarded interaction. In section III we analyze the relax¬ 
ation dynamics of doublons after a photo-doping pulse, the 
doping-induced changes in the electron spectral function and 
in the screening modes, and the effect of the energy distribu¬ 
tion of the carriers on the screening. While most of the results 
pertain to Mott insulators, we also consider the photo-doping 
of strongly correlated metals and in particular the effect of the 
destruction of the quasi-particle peak on the screening. Sec¬ 
tion IV is a conclusion and outlook. 


II. MODEL AND METHOD 
A. U-V Hubbard model 

We consider the single-band U-V Hubbard model on the 
two-dimensional square lattice 

H{t) = -^2 Sij(t)(4*c jc t + h.c.) - n^rii 

(ij)a i 

+ X t/ W( ni T “ DK-I “ |) + '52V(t)(n i - 1 )(nj - 1 ), 

* (»j> 

where Ci a denotes the annihilation operators of a fermion 
with spin a at the lattice site i, rii = + n 4 , s(t) is 
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the hopping integral between neighbouring sites, whose time- 
dependence captures the effect of an in-plane electric field, 
// is the chemical potential, U the on-site interaction energy, 
and V the interaction energy between two electrons on neigh¬ 
bouring sites. The case V = 0 corresponds to the conven¬ 
tional Hubbard model. Using the identity (m — 1 )(ni — 1 ) = 
2(714 —1/2) (rij-f —1/2)+ 1/2 we can combine the interaction 
terms as / ]TV Vijtiifij with n = n— 1 the density fluctuation 
operator (for the case of half-filling), 

Vij — USij + V 5 {ij) , (1) 

and a shift of the chemical potential n —> fl = fx + U/2. 

The grand-canonical partition function is Z = Tr[ 7 ce s ], 
where 7c denotes the countour-ordering operator for the 
Kadanoff-Baym contour which mns from time 0 to some max- 
ium simulation time f max along the real time axis, back to 0 
along the real-time axis, and then to —i /3 along the imaginary- 
time axis. Following Ref.[ 29 ], we express Z as a coherent-state 
path integral, Z = f V[c*, Ci]e s , with the action given by 


and bosonic Green’s functions for this action are 

Gij{t,t') = -\{ci{t)c*(t')) 

with the expectation values defined as (...) = 

± fV[c*,a\[e s ...\. 

B. Nonequilibrium EDMFT 

The EDMFT approximation maps the lattice problem onto 
a single-site effective action 

Sf m -V,C,0] = - i [ dtdt^-^clWGoaMCvit') 

+ + i 0(*Mc(MO™CO}> 

( 6 ) 


S[c*,c] = - ii [ dt^ 2 c* a (t){{-id t - £)<% + s i:j (t)]c ja (t) 

l C ijcr 

+ [• ( 2 ) 

ij ) 

In order to decouple the interaction term we will use the 
Hubbard-Stratonovich identity 



V[x 1 (t),x 2 (t), 


x exp 1 


yj ( 27 r) Ar det 2 l 


( 3 ) 


-'52 x i(t)bi(t)8 c (t,t')} 

i 


where A is a real symmetric positive-definite matrix and b,(t) 
and Xi (t) are fields defined on the contour. We will perform 
the so-called ‘‘U (/-decoupling”,^ where the full interaction 
term is decoupled via an auxiliary bosonic field <^j. Choosing 
bi = i fii, Aij = Vij and Xi = 4 >i the transformed action 
becomes 


S[c*, c, (j)] =—\ J dtdt' 


ijcr 


+ \ fa (*)(*> (t') + i w S c (*> t ') n i{ t ) 

ij i 

( 4 ) 


where we have introduced the fermionic Green’s function 

= [(i d t + fySij - Sij]6 c (t,t'). The fermionic 


whose fermionic (Go(t, t')) and bosonic (ld(t, t')) Weiss fields 
are fixed by a self-consistency condition. This effective action 
is obtained by integrating out all sites but one from the lattice 
action | 4 ]i and taking the infinite dimensional limit. 29 31 The 
hybridization function for the electrons, A a (t, t'), is given by 
Goj(t,t') = [ic? t + ji]Sc(t,t') - A a (t,t') and the equivalent 
bosonic function V corresponds to the retarded component of 
the interaction Id: U(t , t') = U(t)Sc(t, t') + T)(t, t'). In order 
to obtain a purely electronic action we can integrate out the 
bosonic field (f> and obtain: 

jv[<j>}e s ^l=e se e^ u \ ( 7 ) 


where electronic action S e is given by 


* 5 fm P [c*,c] =-i 


dtdt'{ -Y^cl(t)g o j{t,t')c a {t') 

a 


( 8 ) 


In EDMFT, the impurity Dyson equations for Gi mp (f, t') = 
—i (c(t)c*(t')) s e-b and Wj mp = i(</>(f)</>(f , )) s =-6 are given by 


G; mp = Go + Go * £ * Gi m p, (9) 

W iap =U+U*n*W iap , ( 10 ) 

where E (n) is the fermionic (bosonic) self-energy and the 
star denotes the convolution on the contour. The bosonic 
propagator Wi mp and the retarded interaction Id are con¬ 
nected through the charge-charge correlator Ximp = 
(' Tcn{t)n{t')} as (see App. 0> 


H/mp Id Id * Ximp * G. 


(ID 


The solution of the impurity problem, i.e., the calculation 


of Gimp and Wimp is described in Sec. IIC From the Dyson 


equations (| 9 ]» and (|T0ji we obtain the self-energies E, n and the 
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EDMFT approximation identifies these with the lattice self¬ 
energies. The solution of the lattice Dyson equations 

Gk = Go.fe + G 0 ,k * £ * Gk, (12) 

W k =v k +v k *IL*W k , (13) 


then yields an approximation for the lattice Green’s func¬ 
tions G k . W k and from these we can estimate the local lat¬ 
tice Green’s functions G\ oc and Wi oc by averaging over k. The 
EDMFT self-consistency condition demands that G\ oc = Gimp 
and Wi oc = Wi m p. Therefore, updated Weiss fields can be 
obtained from the impurity Dyson equations (|9| and ( fT0| > by 
replacing the impurity Green’s functions with the local lattice 
Green’s functions. The solution of the impurity problem J8 i 
via (11, the impurity and lattice Dyson equations |9]), (QT) i 
and (12 1 , © the EDMFT approximation for the lattice self¬ 
energies and the EDMFT self-consistency equations for G i oc 
and Wioc form the closed set of nonequilibrium EDMFT equa¬ 
tions. 

The nonequilibrium EDMFT calculation is implemented as 
a step-wise time-propagation, which starts from an equilib¬ 
rium EDMFT solution at time t = 0. For each time-step 
along the real-time axis, we iterate the following procedure 
until convergence is reached: 


1. Start with some initial guess for the dynamical mean 
fields A (t,t') and T>(t,t') (for example extrapolations 
of the converged solution for the previous time step). 


2. Solve the impurity problem and obtain Gi mp (f,f') and 


W lrnp (t, t') as described in the Sec. IIC 


3. Obtain a new approximation for the dynamical mean 
fields by closing the lattice self-consistency relations as 
described in the Sec. m 


Since we will use a strong-coupling impurity solver, 17 
the implementation of the self-consistency loop needs to be 
slightly reformulated. While the scheme for the fermionic 
self-consistency loop has been discussed in Ref. [19] we will 
explain the implementation of the bosonic self-consistency 


loop in Sec. IID 


C. Impurity solver 


In order to solve the impurity problem corresponding to the 
action 0 we combine a hybridization expansion and a weak- 
coupling expansion in powers of the retarded density-density 
interaction. It is therefore convenient to first express the parti¬ 
tion function in terms of A and V as Z = l'r, : \Tcc S ), with the 
contour action 

S = - ij J dtdt ' ci(t)A a (t, t')c a (t’) 

+ \j dtdt'n(t)D(t, t')n(t') + J dtH\ oc {t) + const.| 

( 14 ) 


and Uioc(f) = —Jd^2 a n a {i) + U(t)n^(t)ni(t). The double 
expansion then leads to 


OO OO 


2 = EE 

n—0 m—0 

x / dt i... dt. 


(-i)" H)' 


E Tr 


J dt\...dt n ’ J dti ... dt m iTce 1 Sc dtH '<*( t '> 

x c f ai (ii)c CTl (ii). ■ ■ ci n {t„)c an (t' n ) 
x n(fi)n(fi)... n(t m )n(t' m ) 

x A ai • • ■ A CTn (t n , t' n )V(i i, ii)... 


(15) 


In order to evaluate the trace over the electronic configu¬ 
rations one can insert a complete set of eigenstates of H\ oc , 
J2 n \ n )( n \’ between consecutive operators O and factor the 
trace into a product of impurity propagators g and hybridiza¬ 
tion vertices for electrons (F a ) or bosons (K): 

g n (t,t / ) = -i(n\r c e- i f^^\n^ 

Km = (n\c a \m), (16) 

A-nm = ^nm(tr|n|n). 

The Taylor expansion of the partition function can then be rep¬ 
resented as the sum of all diagrams made up of bare impurity 
propagators and vertices connected by A and V lines, in anal¬ 
ogy to the discussion for the Hubbard model in Ref. [19] The 
vertices corresponding to the coupling between the pseudo¬ 
particles and the hybridization function A, and the coupling 
between pseudo-particles and the retarded interaction V are 
shown Fig.|TJa). 

To resum the series we define the pseudo-particle self¬ 
energy E p as a sum of all parts of the above diagrams, that 
cannot be separated into two by cutting one pseudoparticle 
propagator line. With this we can define the renormalized 
pseudoparticle propagator Q via the pseudo-particle Dyson 
equation 


Q = g + 5® ^p®G, (17) 

where the star with arrow denotes the cyclic convolution 
J c dta(t 7 t)b(i,t') restricted to countour-ordered time argu¬ 
ments t > t > t' (cyclic time ordering), for details see Ref. 

E] 

In practice, we need to truncate the self-energy expansion 
at a given order in A and V. To obtain a conserving approx¬ 
imation, we construct a Luttinger-Ward functional $ from all 
vacuum skeleton diagrams involving fully renormalized pseu¬ 
doparticles propagators and the corresponding vertices. The 
lowest order perturbative strong-coupling method is called 
noncrossing approximation (NCA),® 33 because it sums all 
diagrams without crossing A and/or V lines. We will use 
this approximation in the present study for the nonequilibrium 
case. 

A detailed derivation of the strong-coupling equations on 
the Keldysh countour has been presented for the Hubbard 
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a) 


+— = 9 


= A vw = D 


The impurity and the lattice Dyson equations for the bosons 
read 



Figure 1. (a) Vertices representing the coupling between pseudo- 

particle propagators (solid lines) and hybridization functions (dashed 
lines), or between pseudo-particle propagators and retarded interac¬ 
tions (wavy lines), (b) The lowest order diagrams (in A and V ) con¬ 
tributing to the Luttinger-Ward functional, (c) The pseudo-particle 
self-energy diagrams corresponding to the approximation (b) for the 
Luttinger-Ward functional. 


model in Refs. El and [HI In this case, the NCA diagrams 
for *!> contain just one A line. The generalisation to the par¬ 
tial resummation of the series given in Eq. ( p3] > includes an 
additional diagram with a single V line, see Fig. [T] b). The 
expansion of the partition function in Eq. © includes terms 
with crossing electron and boson propagators, but the lowest 
order $ diagram which produces these terms is 0(K 4 F 2 ) or 
0(F 2 K 4 ) and will be neglected in our calculations. 

The self-energies are obtained as a functional derivative of 
the Luttinger-Ward functional with respect to the correspond¬ 
ing pseudoparticle propagators, 

= Wvjy (18) 

and are depicted for the NCA case in Fig. [TJc). The explicit 
expression for these diagrams is 




where the Weiss field is given by U(t, t') = Uo(t, t')+T>(t, t') 
andUo(t,t') = USc (f, t'). Note that the inverse of the bosonic 
Weiss field appears in the Dyson equation, so we have to cal¬ 
culate U~ l . We make the ansatz U~ 1 = Uq 1 — B. Since 
U *U~ 1 = Sc, the components Uq 1 and B have to satisfy the 
condition {Uq + V) * (Uq 1 — B) = Sc- From this and the 
relation Uq 1 = jj5c(t,t') one finds 


5c(t,t)+ jy^V(t,t) 


* B(t, t') 


1 

W) 




1 

W)’ 

( 22 ) 


so that the equation for B(t, t') is a numerically stable integral 
equation. 

Next, we define the bosonic function = 

t') — n(f,f'), which allows us to write the impurity 
Dyson equation as 

W-J e {t,t')=u~ 1 {t, (23) 

(1 + Wi mp * B) *u = W imp . (24) 

In the second line, we thereby obtained a numerically sta¬ 
ble integral equation for u. In solving this integral equa¬ 
tion, instantaneous terms proportional to Sc{t,t') need to be 
treated separately; such terms arise because Wi mp (f,f') = 
U(t)6c(t,t') + W les (t, t') has an instantaneous term. In the 
general case, where the solution may contain a singular term 
we write: [(1 + F s (t))S(t,t) + F re s(t,f)] * [G 5 S c (t,t') + 
G re s(f,f')] = Q s {t)8 c {t,t') + Q res {t,f), where X s (X re s) 
marks the instantaneous (regular) part of the propagator X = 
F. G, Q. Collecting the instantaneous and retarded terms 
yields the two equations 


£ p (f, t') = - i [F^git, t')F a A a (t, £') + F a g(t, t')F a A a (t', t)} G s {t ) = [1 + F 5 (t)]~ 1 Q s {t) 

(19) 


(25) 


D. Lattice Dyson equation 


1 + [1 + F°(f)] - 1 F reg (f, t) 


* G reg (£, t') = 


[1 + F s (t)]- 1 Q reg (t,t') - [1 + F 5 (t)]- 1 F KS {t,t')G 6 (t'). 

(26) 


We need to solve the lattice Dyson equations in order to ob¬ 
tain new approximations for the Weiss functions. In EDMFT, 
the impurity and the lattice Dyson equations for the fermions 
read 

Gfc 1 (£,£') = (i dt + p,)S c (t,t') - E k (t, t') - 

where we have introduced E k [t,t') = e k S{t,t'). In order to 
close the self-consistency loop for the electrons we can pro¬ 
ceed as described in Sec. II.B.4 of Ref. [T9, but the bosonic 
self-consistency loop requires some modifications. 


In the case of Eq. ( |24]> the propagators G, Q and F are G = u, 
Q = Wimp, and F(t{t) = C/(£)B(£,£) + W reg (£,£i)*B(£ 1 ,f), 
which means that F has no instantaneous contribution (F s = 
0). The solution is thus obtained from 

u\t ) = W s (t) = U (t), (27) 

[1 + F}* u reg = WCp(M') - F(t,t')U(t'). (28) 

The lattice Dyson equation can be rewritten 

as W k \t,t') = v^ 1 (t)Sc(t,t') - n(f,f') = 

u+vMtoc (t)Sc(t,t') - n(f,f') = m _1 (£, t’) - A k {t,t’), 

^nonloc 

with A k (t,t') = a k {t)5 c {t,t') and a k (t) = u{u h +v « k (t ) 
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Figure 2. (a) Comparison of the double occupancy in the Holstein-Hubbard model on the Bethe lattice for /3 = 5, U = 10, and u>o = 2 

obtained from the NCA and OCA approximation in combination with the weak-coupling expansion (W) and Lang-Firsov (L) transformation, 
and the numerically exact QMC results, (b) Comparison of the double occupancy for the extended Flubbard model for /3 = 5, U = 12 obtained 
from the NCA approximation and the numerically exact QMC solver on the square lattice. 


This allows us to cast the problem into the form of a stable 
Volterra integral equation for W k (t, t') : 


(1 — it * Ak ) * Wk = u. 


(29) 


We now use again Eqs. ( [25| and ( [26] ) with the substitutions 

G = Wk, Q = u and F — —u * Ak = —U(t)a k {t)5 c (t,t) — 
u reg (f,f)a fc (f): 


W 5 k {t) = 


U(t) 


1 — U (t)ak(t) 


= (U + 


,,nonloc 


l - - 


u les (t,t)a k (t) 


1 — U (t)dk(t) 


m, 

=(M') 


* wr s = - 

1-U(t)a k (t) 

tt res (f,/')a fc ((')(£/ + v k ){t') 
1 — U(t)a k (t) 


(30) 


The sum over k for the instantaneous term W k gives the cor¬ 
rect instantaneous contribution for the local bosonic propaga- 
tor: + v k ° nloc ) = U. In analogy to the electronic case 

we take the sum over k in Eq. (|29]l and use Eq. (t24|) to obtain 
W x \ 


B * Wi„ 


= Y, A k*W k = W U 


(31) 


where in the middle expression the instantaneous contribu- 


nonloc 


) = E 


..nonloc 


k u k 


/u = o. 


tion vanishes, J2k A k(U + v k 
in agrement with the left hand side. After inserting the con¬ 
jugate of Eq. ( [29] ), namely W k = u + W k * A k * u, and 
Wimp = u + Wi mp * B * u into Eq. (31 1 we find 

B + B * Wmp * B = + A k * W k * A k ] = W 2 , (32) 


k 


where the instantaneous contribution to the middle expres¬ 
sion vanishes due to + E k V k U (. U + V k ) = 

J2 k j 7 = 0 . The regular part B of the new U 
be calculated from 


-1 


can now 


[1 + W 1 ]*B = W 2 . 


(33) 


We still need to obtain the expression for the retarded interac¬ 
tion V(t,t') and therefore we once more use Eq. (22 1 in the 
form 


(1 - U(t)B(t, t)) * V(t, t') = U(t)B(t, t')U(t'), (34) 


which is again a stable Volterra integral equation. 


E. Benchmarks in equilibrium 

To test the accuracy of the impurity solver, we compare 
some equilibrium results for the double occupancy with nu¬ 
merically exact Monte Carlo data .- 8 Instead of an EDMFT 
solution for the extended Hubbard model, we first consider 
DMFT results for the Holstein-Hubbard model. In this case, 
there is only one bosonic mode with frequency ■x’q and cou¬ 
pling strength A. The bosonic Weiss held is given by the 
single bosonic mode propagator = A 2 Vo(t,t') = 

—iA 2 Tr[e _1 -Cc dt “ciO +n )/ 2 (j)[t)(f>{t')\, where n denotes the 
conjugate momentum. The black solid and dashed lines in 
panel (a) of Fig. [2] show the double occupation for f3 = 5, 
U = 10, and ojq = 2 obtained from the combined hybridiza¬ 
tion and weak-coupling expansion, with the solid line corre¬ 
sponding to the NCA approximation and the dashed curve 
to the one-crossing approximation (OCA), which considers 
self-energy diagrams with at most one crossing A and/or V 
line. The Monte Carlo result, which considers all relevant di¬ 
agrams, is shown by the green dots. 

For U = 10, the model without phonon coupling is in the 
Mott insulating phase. As A is increased, the Coulomb in¬ 
teraction gets screened and the double occupancy increases. 
Around A ~ 3.0, the solution crosses over to the bipolaronic 
insulating phase, which is marked by a large value of the dou¬ 
ble occupancy, see inset of Fig. [2]a). (At lower temperature, 
the transition to the bipolaronic insulator occurs via an inter¬ 
mediate metallic phase, see Ref. (28]) As expected, the NCA 
approximation overestimates the correlation effects and thus 
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Figure 3. (a) Equilibrium spectral functions for f3 = 20 and indicated values of U and V. (b) Equilibrium phase diagram in the space of U 
and V at/3 = 20. The dashed line roughly indicates the metal-insulator crossover defined by the appearance of a quasi-particle peak, (c)-(d) 
The equilibrium spectral function A(lj) for [7 = 12 (c) and U = 7 (d) for indicated values of the nearest neighbor interaction V. The insets 
show the spectra on a logarithmic scale. 


underestimates the double occupation, while the OCA approx¬ 
imation is quantitatively more accurate. We also see that the 
weak-coupling treatment of the electron-phonon interaction 
correctly captures the screening effect, i.e. while the double 
occupancy is shifted to lower values (as expected in an NCA / 
OCA calculation), it exhibits the correct A-dependence in the 
weak-coupling regime. The method however does not capture 
the transition into the bipolaronic state. 

For comparison, we also show with red lines the equilib¬ 
rium results from an alternative scheme in which the phonons 
are first decoupled by a Lang-Firsov transformation and then 
integrated out. The Monte Carlo method is actually based on 
such a decoupling, 27 ^ and the NCA/OCA approximation of 
this Lang-Firsov approach has been discussed in Ref. [34] In 
the small-A regime, this approximate method is slighlty less 
accurate than the combined hybridization and weak-coupling 
expansion discussed in this work. On the other hand, the 
Lang-Firsov method is not limited to small phonon couplings 
and correctly captures the crossover to the bipolaronic phase. 

We next test the EDMFT solution for the U-V Hubbard 
model. The comparison of the double occupancy rid to the 
QMC results shows a similar trend as found in the Holstein- 
Hubbard calculations, see Fig.[2]b). The NCA approximation 
underestimates the double occupation, but correctly captures 
the screening effects in the weak coupling regime V < 3.5. 


Since the combined NCA and weak-coupling approach pro¬ 
vides a qualitatively correct description of the physics in the 
Mott insulator, as long as one stays away from the charge or¬ 
dered phase, we will use it in the rest of this work to inves¬ 
tigate the real-time dynamics of photo-doped Mott insulators 
and of models in the metal-insulator crossover regime. 

III. RESULTS 

We use the nonequilibrium EDMFT scheme to simulate the 
time evolution of a U-V Hubbard model in which an elec¬ 
tric field pulse produces a nonthermal occupation of doublons 
and holes. The photo-doping induced changes and the relax¬ 
ation and eventual thermalization are illustrated by measur¬ 
ing the fermionic and bosonic spectral functions, as well as 
the double occupancy. To incorporate the electric field into 
the model, we use a gauge with pure vector potential (A), 
so that the electric field is given by E = —<9 t A. The vec¬ 
tor potential enters Eq. 0 via the Peierls substitution, i.e. 
the hopping integrals Sij(t) acquire a time-dependent phase 
factor, or, eq uivalently, the band energies et are shifted as 
(k Specifically, we use a pump pulse of the form 

E(t) = E 0 sm(uj(t — f 0 )) exp(— 4.6(f — fo) 2 Ao) with f re " 
quency u, amplitude Eq : and a Gaussian envelope. The width 
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Figure 4. (a),(b) Relaxation dynamics of the double occupancy for different values of V after a pulse with frequency w = 10 and pulse 
amplitude Eo = 5. The on-site interaction is U = 10 (a) and U = 7 (b). 


of the pulse, to = 2tt/u>, is chosen such that the envelope ac¬ 
commodates a single cycle. Unless otherwise stated the pulse 
frequency is oj = 8.0. The calculations are done for a 2D 
square lattice with bandwidth W = 8|s| and the field is point¬ 
ing along the lattice diagonal. We use |s| as the unit of energy, 
and h/\s\ as the unit of time. 


A. NCA phase diagram 

Before proceeding to the nonequilibrium dynamics we will 
consider the equilibrium properties in two different phases, 
namely the paramagnetic Mott insulating and metallic phase. 
In equilibrium the Mott insulator is characterized by a well 
defined gap in the equilibrium spectral function A(u>) = 
-ilm[G fl (w)], while the metallic phase shows a coherent 
quasi-particle peak, see Fig. [3ja). Within the NCA approxi¬ 
mation we are limited to rather high temperatures, T > 1 /20. 
Although this temperature lies above the end-point of the 
Mott transition line, one nevertheless still observes a relatively 
sharp crossover between the metallic and insulating regimes. 
We determined the boundary between the metallic and insulat¬ 
ing regimes by the (dis)appearance of the quasi-particle peak 
in the equilibrium spectral function. This crossover line is 
plotted in Fig. [3jb). 

In order to understand the effect of the nonlocal interac¬ 
tion V we show the spectral functions for U = 12, deep in the 
Mott insulator, in Fig.[3jc) and for a system close to the MIT at 
17 = 7 in Fig.[3|d). In the Mott insulator, increasing the near¬ 
est neighbor coupling V causes a slight shift of the Hubbard 
band toward lower frequencies and an enhanced weight of the 
high energy plasmon satellite, see inset of Fig. [3jc). Note that 
the high-energy peak is already present in the V = 0 case 
(Hubbard model), where it is a consequence of higher order 
processes entering the strong-coupling diagrams via the hy¬ 
bridization function. 37 This can be easily understood in the 
case of the Bethe lattice, where the lattice self-consistency 
condition simplifies to A = |s(0)| 2 G; mp (f, t'). The 
hybridization function is within NCA given by a bubble of 
pseudo-particle Green’s functions,® which also leads to ex- 


Table I. Values of the reduced effective interaction U{uj = 0) mea¬ 
sured at t = 20 after a pulse of frequency w = 10 and amplitude 


= 5. 

V 

U = 10 

U = 7 

0.5 

9.97 

6.92 

1.0 

9.86 

6.67 

2.0 

9.29 

5.47 


citations of 2U from the lower Hubbard band, i.e. to a side¬ 
band at uj ss —U/2 + 2U = 317/2. The enhancement of this 
peak with increasing nearest neighbor interaction V is a con¬ 
sequence of plasmonic excitations, with characteristic energy 
17, from the upper Hubbard band. These excitations contribute 
spectral weight at the same energy w ss 17/2 + U = 317/2. 
Close to the MIT the increased screening due to the inter-site 
interaction leads to a crossover to a coherent metallic state 
as a function of V (cf. data for V = 2 and 17 = 7 in 
Fig.Rld)), and a significant increase of spectral weight at high 
freq^ncies.™! 

B. Enhanced doublon relaxation due to dynamical screening 

In the following we focus on the nonequilibrium dynamics 
deep in the Mott insulating regime and close to the crossover 
line. In Fig. [4] we plot the time evolution of the double occu¬ 
pancy after the pump. The amplitude of the pump is chosen 
such that a relatively high density of charge excitations is cre¬ 
ated in the system. (For 17 = 10, V = 2 the change in the 
double occupancy is A rid ~ 0.04. All the results shown in 
Fig .|4]are computed for a fixed amplitude Eq, and thus corre¬ 
spond to similar values of An,j, because the double occupancy 
only weakly depends on the nearest neighbor interaction V.) 

After a transient dynamics the double occupancy follows 
an exponential relaxation. In the Hubbard model, the relax¬ 
ation time increases exponentially with increasing U if the 
gap size is large.®® In small gap insulators, impact ioniza¬ 
tion processes can lead to a rapid carrier multiplication at 
short times, followed by a slower exponential relaxation.® As 
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Figure 5. Relaxation time of the double occupancy as a function of nearest neighbour coupling V for U = 10 (a), U = 7 (c), pulse frequency 
lj — 10 and pulse amplitudes Eq = 2.5 and 5. The dashed lines represent the relaxation time in the Hubbard model with instantaneous 
interaction U = U(ui = 0) corresponding to the static interaction in a U-V Hubbard model with a given V (see discussion in the main text). 
(c,d) Scaling of the difference between the inverse relaxation times at V = 0 and V > 0 with V 2 for the same parameters as in panels (a,c). 


can be seen from the equilibrium spectra plotted in Fig. [3ja), 
even U = 10 lies within this small gap regime, so that we 
will focus on the long-time relaxation. The relaxation dy¬ 
namics for different values of the nearest neighbour inter¬ 
action V are plotted in panels (a) and (b) for U = 7 and 
U = 10, respectively. We fit the relaxation by an exponen¬ 
tial function ndit) = rid{t = oo) + Aexp(—f/r) in the range 
5 < t < 30 to extract the relaxation times. While the relax¬ 
ation curves look qualitatively similar, the relaxation times, 
plotted in Fig. [5ja) and (c), are strongly reduced for larger V. 


Several effects could potentially explain this observation. 
On the one hand, the dynamical screening reduces the effec¬ 
tive instantaneous interaction U(u = 0) (see Sec. HID for 
a detailed analysis), so the dynamics may be more properly 
described by a Hubbard model with an effectively reduced re¬ 
pulsion U = U(oj = 0), which in turn leads to a faster relax¬ 
ation. In order to test this scenario we measured the screened 
effective interaction U(uj = 0) after the photo-doping pulse 
for different values of the nearest neighbor interaction V, see 
Table [IJ and then recalculated the dynamics in the Hubbard 
model (V = 0) using these reduced interaction parameters. 
The corresponding relaxation times are shown by dashed lines 
in Fig.[5ja) and (c). At U = 7, in the metal-insulator crossover 
regime, the reduction in the effective on-site interaction (and 
hence gap-size) contributes substantially to the faster relax¬ 


ation times, but in the Mott insulator case (U = 10), this effect 
cannot explain the observed large changes. 

The second scenario is that the relaxation rate is enhanced 
due to the absorbtion of bosonic collective excitations (plas- 
mons), which opens an additional relaxation channel. The 
coupling to the bosonic excitations can be understood from 
the form of the electron-boson action in Eq. (|6j, which is 
equivalent to the action of a Anderson-Holstein model with 
a coupling to a continuum of free bosons. 29 Previous studies 
of the relaxation dynamics in the Hubbard-Holstein model 12 ® 
showed an enhanced relaxation of the charge excitations due 
to the coupling with phonons. In order to check the lat¬ 
ter idea we can use the empirical Matthiessen’s rule 1/r = 
1/tv=o + l/i~b to separate the inverse relaxation time into an 
electronic and bosonic contribution. Based on this formula we 
extract the inverse relaxation time due to the scattering with 
the bosonic bath 1 /tj, from the relaxation times obtained by 
fitting. The result is presented in Fig. [5jb,d). 

We note that the extracted relaxation times 1 /rj are propor¬ 
tional to V 2 for small V. One can understand this scaling 
from the effective Hamiltonian representation of the impu¬ 
rity model, which corresponds to an Anderson impurity model 
coupled to the bath of bosonic modes. 29 ^The relaxation time 
due to the coupling to the bosonic bath can be approximated 
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Figure 6. (a,b) Partial Fourier transform of the spectral function A(w, t ) with insets showing the difference from the initial value A(uj, t) — 
A(lo, t = 0) for U = 10, V = 2,Eo = 5.0 (panel (a)) and U = 7, V = 2, Eq = 3.5 (panel (b)). The dashed black lines represent the initial 
equilibrium spectral function. (c,d) Change in the gap size measured by the frequency lo* at which the spectral function is equal to some fixed 
value, A(uj*) = Aq, for U = 10 (c) and U = 7 (d). The dashed lines represent the u>* for the same parameters but with coupling to an 
external bath with A = 1, see also the discussion in the main text. For U = 7, there is no intersection for Ao = 0.05 in the case of coupling to 
external bath. The corresponding cuts through the spectral function are shown by the dashed lines in panels (a,b). 


by the Fermi golden rule: 


C. Spectral properties 


y - Uk) ~ A(lu f ) ^ \ 2 k — 

A( k (35) 

— du\m\D(uj)\ oc A(u}f)V 2 , 


where A;,,, LOk are the coupling constants and phonon frequen¬ 
cies for the fc-th bosonic mode and A(u>) is the electronic 
spectral function. The main simplification used was to ap¬ 
proximate the spectral function at the final energy A{u>f = 
(jj — ujt-) as a constant. The numerical investigation in Ref. [29 
showed that for weak nearest neighbour interaction V the in¬ 
tegral over the bosonic Weiss field scales as f dwlm[2?(u;)] a 
V 2 . Hence, our phenomenological analysis shows that the en¬ 
hancement of the relaxation due to the nonlocal interaction V 
and the corresponding V 2 scaling is consistent with a relax¬ 
ation aided by the coupling to a bath of bosonic degrees of 
freedom (plasmons).EE3El 


Further insight into the relaxation dynamics is obtained 
from the evolution of the partially Fourier transformed spec¬ 
tral function A(t, to) = — ^Im dt', t), 

where we use f max = 10. First we study the effect of the 
screening in the Mott insulating phase (U = 10), where the 
equilibrium spectral function consists of the upper and lower 
Hubbard band separated by a well defined gap (Fig. [6jr). The 
pulse leads to an increase in the number of charge carriers 
(doublons and holes), which in turn results in stronger screen¬ 
ing. The spectral weight at cu = 0 increases slightly due to the 
effect of screening, heating and doping. The enhanced screen¬ 
ing effect manifests itself on the timescale of 1/bandwidth 
(note the nonuniform time mesh in Fig. [6ji. At longer times 
the screening is further increased due to the redestribution of 
the charge carriers within the upper Hubbard band and on even 
longer time scales due to doublon production associated with 
thermalization. 

In order to compare these results to the evolution of a sys¬ 
tem close to the metal-insulator crossover (U = 7), we adjust 
the pulse strength E 0 such that the same amount of photo- 
doped charge carriers 5rid « 0.04 is present after the excita- 
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tion. The result is shown in Fig.[6jb). While the quasi-particle 
peak quickly disappears there is again a shift and broadening 
of the Hubbard bands, associated with an increase in spectral 
weight in the pseudo-gap region. The insets in Fig. [6ja),(b) 
show the evolution of the difference of the time-dependent 
spectral function A(ui,t) — A(ui,t = 0) for w > 0. The 
spectral weight at the lower and upper edge of the Hubbard 
band is increased and as a consequence the gap is reduced and 
partially filled in. Since the system is almost thermalized for 
t > 10 these changes can be attributed to the heating of the 
system, which we have confirmed by comparison to the equi¬ 
librium spectral functions at elevated temperatures. 

In particular, the increase of weight at low w is related to 
the partial filling-in of the gap, while the increase at large u is 
related to the appearance of side-bands in the photo-doped or 
thermally excited system (See Fig. [3jc, d)). 

In order to analyze the doping and screening induced 
changes in the gap size we plot the evolution of the fre¬ 
quency w*, where the spectral function takes some fixed 
value: A(u>*) = A 0 , see Fig.[6]p,d). The corresponding cuts 
are shown as dashed lines in Fig. [6ja,b). At the shortest times 
the charge carriers are inserted near the middle of the upper 
Hubbard band, which results in a distortion of the spectral 
function (the dashed black line shows the initial equilibrium 
result). While the pulse is on (t < 3) there are field induced 
effects, which lead to a nontrivial dynamics of the spectral 
function. The doping, heating and screening induced changes 
result in an asymmetric reshaping of the Hubbard bands (see 
insets) and the modifications in the gap size are quite dif¬ 
ferent from what one would obtain from a rigid shift of the 
Hubbard bands of the undoped system. In order to investi¬ 
gate the reduction of the gap quantitatively we analyze the 
difference of the frequency w* from the equilibrium value. 
Aw* (t) = to* ( t ) — oj* q . The rapid reduction of Aw* just after 
the excitation is related to the fast reduction of the effective 
U due to screening, see Sec. Ill D The longer time dynam¬ 
ics depends on the gap size: for U = 10 the gap keeps de¬ 
creasing even on the longest accessible times scales, while for 
U = 7 the gap size rapidly stabilizes. This difference can be 
explained by the longer thermalization time in the C7 = 10 
case. 

The dashed lines in Fig. |6fc,d) show the analogous results 
for a model with local coupling to a heat bath; for details and 
further analysis see Sec |IV| Due to the cooling by the heat 
bath, the gap recovers its original value, as will be further dis¬ 
cussed in Sec. IIV Cl 


D. Effective interaction 

Within EDMFT, the inter-site interaction translates into a 
retardation of the effective on-site interaction U(t,t'), while 
the fully screened interaction W (t. t') also includes the lo¬ 
cal screening effects. In equilibrium and in the Mott insulat¬ 
ing phase the imaginary part of the fully screened interaction 
Im[W(w)] consist of one broad peak around w « U, which 
is weakly shifted to lower frequencies upon increasing y.USSO 
In the chemically doped Mott insulator, or in the strongly cor¬ 


related metal phase, Im[FF(w)] exhibits a two-peak structure. 
The low energy peak is associated with transitions between 
quasi-particle peak and Hubbard bands, while the higher en¬ 
ergy peak at w ~ U is associated with transitions between the 
Hubbard bands. By increasing the doping, the weight of the 
low energy peaks is strongly enhanced. The real part of the 
fully screened interaction Re[W(w)] is reduced for low fre¬ 
quencies as we increase the nearest neighbour interaction V 
or the doping, reflecting the stronger screening effect. 

An interesting question is how these spectral functions 
change after a photo-doping excitation in a Mott insula¬ 
tor. In order to get insights into the nonequilibrium dy¬ 
namics of screening we perform the partial Fourier trans¬ 
form of the fully screened interaction to obtain W R (t,u)) = 
j-t+W dt > e -Mf-t) W R(f^ m T he results for U = 10, V = 
2.0 and a photo-doping concentration of A rid = 0.04 are 
shown in Fig. |7](a),(c),(e). In equilibrium the imaginary part 
of W shows a broad peak at w p ~ U, while the real part 
is slightly reduced for w < w p as a result of local screen¬ 
ing. During and shortly after the pulse (t < 3) there is a 
strong modification in the distribution and coupling strength 
of the screening modes. Due to the photo-doping and the 
screening from low-energy excitations within the Hubbard 
bands, a broad low-energy peak appears in the imaginary 
part of W R (t, ui) and the instantaneous effective interaction 
Re[VF i? (f, oj = 0)] is strongly reduced. 

Even though the screened interaction W{ui) of the photo- 
doped Mott insulator looks qualitatively similar to that of a 
chemically doped Mott insulator, the origin of the low-energy 
peak in Im[FF fl (w)] is different: The single-particle spec¬ 
tral function of photo-doped Mott insulators does not feature 
a quasi-particle peak at oj = 0, nor sharply defined quasi- 
particle features near the gap edges.® The low-energy peak 
in Im \W R (u))\ is therefore not associated with transitions be¬ 
tween a well-defined quasi-particle band and the Hubbard 
bands, but rather with excitations within the Hubbard bands. 
As a result, the feature is broader in the photo-doped system 
than in a chemically doped one. 

A useful quantity to characterize the screening is the aver¬ 
age boson frequency 


J^duujlm[W{t,uj)} 

/ 0 °° duj ImflE (t , w)] 

It strongly decreases during the pulse (t < 2.5), as shown in 
the inset of Fig. [TJa.b). In the case of U = 10 the initial fast 
drop of TD{t), which is a consequence of the doping-induced 
appearance of the low energy mode, is followed by a slower 
long time relaxation associated with changes in the energy 
distribution of the photo-carriers. For (7 = 7, the system is 
essentially thermalized after t ~ 7 and no further changes in 
the bosonic spectral function or average screening frequency 
occur. 

We note that the changes in the screened interaction are less 
dramatic for U = 7 than for U = 10, see Fig. [7ja,b). This 
can be understood by the fact that the initial state is metal¬ 
lic, and already has low-energy screening modes. The heating 
destroys the quasi-particle peak (i.e the system moves across 
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Figure 7. Real (a,b) and imaginary (c,d) part of the partial Fourier transform of the screened interaction Re[VF(f, u>)\ — U for U = 10, Eo = 5 
(left column) and U = 7, Eo = 3.0 (right column) at fixed V = 2,lo = 8. The insets in (a,b) show the evolution of the average boson 
frequency w(f). The black dashed lines represent the screened interaction in the initial equilibrium state. (e,f) Real part of the partial Fourier 
transform of the regular part of the partially screened on-site interaction Re[27 rcg (t, w)] = Re[W(f, w)] — U for U = 10, Eo = 5 (e) and 
U = 7, Eo = 3.5 (f). 


the metal-insulator line in temperature), so that the final state 
is a thermally excited insulator for which the screening is not 
much larger than in the metallic initial state. In contrast, for 
U = 10 the initial and final states correspond to a cold and 
hot Mott insulator with a very different number of thermally 
excited carriers. The partial Fourier transform of the regular 
part of the effective on-site interaction D(t,t') = U(t,t') — 
1 5 c (t,t')U shows a behavior which is qualitatively similar to 
that of the screened interaction W(t. t'). The main difference 
is that the reduction is smaller (see Fig. [7|e,f)), since the lo¬ 
cal screening effects are absent. The time-dependent changes 
in the effective interactions are consistent with the observed 


changes in the spectral function, where a strong reduction of 
the gap size is observed on the time scale 1 /bandwidth. 

Photo-doping a metallic state destroys the quasi-particle 
peak, 39 which leads to a loss of low-energy excitations and 
hence low-energy screening. This effect competes with the in¬ 
crease of the screening due to the photo-doping. In a strongly 
excited system, the two effects cannot be easily disentangled. 
To single out the effect of the destruction of the quasi-particle 
peak, we choose a low-frequency and low amplitude pulse, 
namely o>o = 2.0 and Eq = 0.25, which reduces the photo¬ 
doping effect. With this pulse we indeed observe a decrease of 
the screening in the low energy regime w < 2 of the real part 
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Figure 8. Real (a) and the imaginary (b) part of the partial Fourier transform of the screened interaction W(t,u) for U = 6.9, V = 2, 
E = 0.25, and u> = 2. (c) Comparison of the dynamics of the partial Fourier transform of the spectral function A(oj, t) for the U-V (full lines) 
model and the Hubbard model (dashed lines) with a correspondingly reduced on-site interaction U = 6.425. 


of the partial Fourier transform of the fully screened interac¬ 
tion W ( t , w), see Fig. [8] In addition to this increase of the real 
part of W, the low-energy feature associated with excitations 
within the quasi-particle band disappears. The latter excita¬ 
tions are responsible for the pole-like structure near io = 0.3, 
which disappears rapidly in agreement with the dynamics of 
the spectral function, see Fig. [8jc) and inset of (a). 

The reduced screening effect feeds back onto the spectral 
function, and leads to a further decrease of spectral weight in 
the low-energy region in comparison with the dynamics of the 
Hubbard model with appropriately reduced U. For the com¬ 
parison in panel (c), we have chosen U such that the equilib¬ 
rium spectral function reproduces the result of the U-V Hub¬ 
bard model at low-energies (u> < 2) in the initial state. As 
can be seen, the reduction of the spectral weight is more pro¬ 
nounced in the case of the U-V model due to the change in 
the screening environment. 


IV. COUPLING TO A THERMAL BATH 

To describe the dissipation of energy to external degrees of 
freedom we need to couple the system to some environment. 19 
On short times, the dissipative environment will lead to a re¬ 
distribution of the energy of the photo-doped carriers (intra¬ 
band relaxation), while on longer timescales the bath en¬ 
hances recombination processes, and enables the system to 
relax back towards the initial equilibrium state. Here we will 
in particular study how the screening is influenced by the 
transient modification of the energy distribution of the photo- 
doped carriers. 

Technically, one can integrate out the environment to ob¬ 
tain an effective description of the system by adding a corre¬ 
sponding contribution to the electronic self-energy E. We use 
the lowest order diagram for a Holstein like electron-phonon 
coupling,® namely Edj ss (w) = A 2 G(t, t'), where A is 

the coupling strength and Bo(t, t') is the free bosonic propa¬ 
gator with a linear spectral density between a low-energy and 
high-energy cutoff, namely — Alm[i?o(u;)] = u>, if uii < oj < 


ujh with u>i = 0.2 and uih = 1.0, and zero elsewhere. The 
bath temperature is set to T = 1/20. Since this treatment is 
only suitable in the weak coupling regime we will restrict our¬ 
selves to values of the coupling (A < 1.0), for which the bath 
only results in a slight broadening of the spectral function in 
equilibrium. The spectral distribution of the bath has a sub¬ 
stantial effect on the relaxation dynamics. To efficiently re¬ 
lax the states around the quasi-particle peak the system needs 
to be able to excite low energy bosons, while a restriction to 
only low energy bosons will slow down the relaxation of high- 
energy carriers due to the necessity for high order processes. 
For this reason, we choose a continuous spectrum which in¬ 
cludes low as well as relatively high energy bosons. 


A. Relaxation with bath 

In contrast to the system without bath, which is approach¬ 
ing a thermal equilibrium state at higher temperatures, the ad¬ 
dition of the heat bath ensures that the system relaxes back to 
the initial state, as can be seen in the evolution of the double 
occupancy in Fig.[9ja,b). This behavior is consistent with the 
previous investigation of the Hubbard model. 39 The question 
which we would like to address here is the effect of the bath 
on the screening. 

For the Mott insulator case (7 = 10 the bath tends to reduce 
the screened interaction Re[TF(w)] at low frequencies on the 
observed time scale. This is a consequence of the faster re¬ 
laxation of the highly excited doublons to the lower edge of 
the Hubbard band, which enhances the screening, compare 
Figs. [9jc) and[7ja). The increase in Re IT'(a,') at later times is 
the result of the bath-enhanced doublon-hole recombination. 
Deeper in the Mott phase (U = 14), the recombination is sup¬ 
pressed and even with a coupling to a heat bath the double 
occupation remains approximately constant on the accessible 
time scales (not show). The bath nevertheless enhances the re¬ 
laxation of the high energy doublons to the lower edge of the 
upper Hubbard band, which leads to monotonically increas¬ 
ing screening effect with time, see Fig.^e). For (7 = 7, the 
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Figure 9. Relaxation dynamics of the double occupancy rid for U = 7 (a), U = 10 (b), and U = 14 (c) with different values of the coupling 
to the bath A during and after the pulse with frequency u = 8 and pulse amplitude Eo = 5. The real part of the partial Fourier transform of 
the screened interaction Re[VK reg (f, cu)] = Re[W (t, w)] — U for U = 7 (d), U = 10 (e) and U = 14 (f) at fixed V = 2.0, A = 1.0, Eo = 5. 


screening dynamics is non-monotonous: the initial increase 
of charge carriers enhances the screening, but the coupling to 
the bath leads to a faster recombination of the charge carriers, 
which eventually reduces W(u,t) to a function close to the 
initial equilibrium screened interaction, see Fig. [9jd). 

B. Photo-doping in the metal-insulator crossover regime 

An appealing scenario would be the appearance of a tran¬ 
sient metallic state (with quasi-particle peak) after the photo¬ 
doping of an insulating initial state, as a result of enhanced 
screening. To investigate this possibility, we have systemati¬ 
cally studied the photo-induced dynamics for a system close 
to the metal-insulator crossover line, where this effect may be 
expected to occur. 

We found that photo-doping an insulator in the vicinity of 
the metal-insulator crossover line can trigger a nontrivial evo¬ 
lution which reflects the effects of increased scattering and of 
the changes in the screening environment, but we were not 
able to realize a screening-induced metal state even close to 
the metal-insulator regime, where it is easy to realize a set-up 
in which U{u> = 0) drops below the critical static value for 
which one finds a metallic solution in equilibrium. A poten¬ 
tial reason may be the detrimental effect of doublon/hole scat¬ 
tering on the emergence of coherent quasi-particles. This re¬ 


sembles the observation that the build-up of a coherent quasi¬ 
particle peak in a photo-doped Mott-insulator takes at least 
longer than the quasi-particle lifetime. 39 We also used rather 
weak pulses, since stronger pulses enhance the destruction of 
the quasi-particle peak and in addition lead to stronger heat¬ 
ing. Another limitation is the NCA based impurity solver, 
which systematically underestimates the metallic nature of 
the solution, and which does not allow us to study the low- 
temperature behavior in the vicinity of the first-order metal- 
insulator transition. 

Photo-doping in the vicinity of the metal-insulator 
crossover line can nevertheless trigger a nontrivial evolution 
which reflects the effects of increased scattering and of the 
changes in the screening environment. Figure [10] shows the 
change in the kinetic energy and in the double occupation af¬ 
ter a weak pulse in systems with weak and strong coupling to a 
heat-bath. We compare the system at U = 7.1, where in equi¬ 
librium a weak quasi-particle peak is present and at U = 7.3, 
where the system only exhibits a pseudogap. Let us first con¬ 
sider the insulating system with weak coupling to a heat-bath 
(A = 0.5). After the pulse, the kinetic energy drops, while 
the doublon number increases, as expected in a photo-doped 
insulator with dissipation. At longer times, the recombination 
of doublon-hole pairs leads to a slow relaxation back to the 
initial equilibrium state. The system with strong coupling to 
a heat bath (A = 1) shows a more puzzling behavior. After 
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an initial drop, the kinetic energy starts to increase beyond the 
initial equilibrium value, while the double occupancy drops 
below the value in the initial state. Both observations suggest 
that a transient state is induced which is even more insulating 
than the initial state, as a result of reduced screening and an 
enhanced effective interaction. 

This is consistent with the dynamics of Re[Z?(w)] which is 
increased at low energies for at all times, see Fig. |T(j|e), so 
that electrons are moving in an effectively more insulating- 
like system. Similarly the spectral function at w = 0 shows an 
initial decrease, which is followed by a slow re-filling of the 
pseudo-gap (see Fig.[I()lc,g)). 

The opposite procedure is to start in an initially metal¬ 
lic state and by applying an electric field pulse destroy the 
quasi-particle peak, which results in a reduced screening ef¬ 
fect. Photo-doping an initially metallic state also leads to an 
increase in AE^o (reduction in the absolute value of the ki¬ 
netic energy), and after a brief transient a drop in the number 
of doubly occupied sites. The time resolved spectral function 
shown in panel (d) reveals a rapid destruction of the quasi¬ 
particle peak, which is followed on longer time-scales by a 
slow recovery. The destruction of the quasi-particle peak re¬ 
duces the low-energy screening, see Fig.fTOff), since conduct¬ 
ing electrons are removed from the system. Despite the slow 
reappearance of the quasi-particle peak, the values of E^ n and 
rid a t the longest times indicate that the system is consider¬ 
ably more strongly correlated than in the initial state, so that 
the system may again be regarded as trapped in a transient in¬ 
sulating nonequilibrium state. Despite the strong coupling to 
the heat-bath, the relaxation back to the inital metallic state 
appears to be inhibited on the accessible time-scales, presum¬ 
ably due to strong doublon-hole scattering. Note that this ef¬ 
fect was already present in the system studied in the previous 
subsection, where the double occupancy rid for U = 7, A = 1 
was reduced below the initial value, see Fig. |9jb), but there, 
due to the competing effect of the stronger photo-doping, the 
dynamics of screening was non-monotonous. In order to elim¬ 
inate this competing effect we used here a weak pulse, which 
prevents a substantial photo-doping, but still results in the de¬ 
struction of the quasi-particle peak. 


lution of the spectral function A(t, ui) for the Hubbard model 
(V = 0) and the U-V Hubbard model with V = 2. The den¬ 
sity of photo-doped carriers is A rid ~ 0.04. At the shortest 
times the charge carriers are inserted near the middle of the 
upper Hubbard band. While the pulse is still on, there is a 
strong distortion of the spectral function due to field induced 
effects. The inclusion of the thermal bath enhance the sub¬ 
sequent relaxation of the high energy doublons to the lower 
edge of the Hubbard band. This leads to the formation of a 
pronounced shoulder-like feature at the band edge both for 

V = 0 and V = 2, in agreement with previous results for the 
Hubbard model. 39 In contrast to the Hubbard model, however, 
which shows only a formation of the shoulder-like feature, for 

V = 2 the increased screening effect leads to a reduction of 
the gap, see Fig. |TTfa). 

Following the analysis in Sec. IIIC we analyze the fre¬ 
quency w*(i), where the spectral function takes some fixed 
value: A(i,w*(f)) = Ay. The reduction of the gap Aw* is 
defined as the difference between the long time average w* 
(measured between times t = 15 and t = 25) of the excited 
system and the equilibrium system Aw* = uj*(Eq) — w* q . 
The reduction of the gap Aw* versus the change in the double 
occupancy Arid for different pulse strengths between 4.0 < 
Eq < 15.0 is shown in Fig. 1 1 1 \c). The model with screening 

V = 2 shows a decrease of the gap for all cuts at different 
energies. In the case of low photo-doping A rid < 0.015 the 
decrease of the gap is linear. The fast decrease and the later 
saturation for Aq = 0.04, 0.05 is a consequence of the for¬ 
mation of the shoulder like feature at the edge of the Hubbard 
band, while at lower energies A = 0.01, w* monotonously 
decreases with increasing density of the photo-doped carri¬ 
ers, due to the partial filling-in of the gap. Therefore in all 
of the analyzed cases the effect of screening is stronger than 
the distortion of the spectral function due to photo-doping. 
The Hubbard model shows a quite different behavior: the ini¬ 
tial reshaping of the Hubbard band leads to an increase of 
the Hubbard gap for low energy cuts ( A 0 = 0.01), while for 
Aq = 0.04,0.05 the initial decrease in Aw* is a consequence 
of the formation of the shoulder like feature. The increase 
at even higher photo-doping is a consequence of the stronger 
reshaping of the upper Hubbard band. 


C. Reduction of the gap 

The effect of the screening on the reduction of the gap in 
the presence of the thermal bath was presented in Fig. [ 6 jc)- 
(d) by the dashed lines, which indicate a gradual recovery of 
the gap back to its original size. In the case of U = 10, see 
Figjhjc), and the cut at A 0 = 0.01 we see the decrease and 
subsequent reappearance of the gap, while the cut at higher 
energy (A {) = 0.03) shows a decrease on the longest availble 
time scales. For the metallic solution U = 7 the gap is also 
initially reduced, but for longer times it starts to recover. The 
different behavior of the two cases is the consequence of the 
longer thermalization times in the insulating case U = 10. In 
order to eliminate the effect of charge recombination we will 
study a system deep in the Mott insulator, namely U = 14. 
To present the effect of screening we will compare the evo- 


V. CONCLUSIONS 

We have discussed the nonequilibrium generalization of ex¬ 
tended dynamical mean field theory, which allows to study the 
screening effect in correlated lattice systems in real time. To 
solve the nonequilibrium EDMFT equations, we have intro¬ 
duced a perturbative impurity solver which combines a self- 
consistenstent hybridization expansion with a weak-coupling 
expansion in the retarded interaction. This method should 
yield qualitatively correct results in Mott insulators or the 
metal-insulator crossover regime at elevated temperature, and 
for systems where the effect of screening is relatively weak 
(no tendency of charge-ordering). 

The formalism has been applied to the half-filled U-V Hub¬ 
bard model on the square lattice, which is driven out of a 
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Figure 10. Relaxation dynamics of the kinetic energy E^m (a) and the double occupancy rid (b) for U = 7.3 and U = 7.1 at V = 2.0 for 
different values of the coupling with the bath A = 0.5, 1.0 during and after a pulse with uj = 8 and Eo = 1.0. The partial Fourier transform 
of the spectral function A{uj, t ) for U = 7.3 and U = 7.1 is shown in panels (c) and (d). The dashed lines corresponds to the equilibrium 
spectral function. Panels (e,f) plot the real part of the partial Fourier transform of the regular part of the partially screened on-site interaction 
Re[75(t, w)] = R e[U(t, w)] — U for U = 7.3 (e) and U = 7.1 (f). The time evolution of the spectral weight A(uj = 0, f) for U = 7.3 (g) and 
U = 7.1(h) for A = 1. 
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Figure 11. (a,b) Partial Fourier transform of the spectral function A(uj,t) for U = 14, V = 2, Eo = 10.0, wo = 14 (a) and U = 14, V = 
0, Eo = 10.0, w = 14 (b) and coupling to a bath with A = 1.0. The dashed black lines represent the initial equilibrium spectral function. 
Change in the gap size oj* (see discussion in the main text) as a function of the change in the double occupancy A rid for V = 2 (c) and V = 0 
(d). 


strongly correlated equilibrium state by a single-cycle laser 
pulse. The photo-doping of carriers into a Mott insulator 
leads to a rapid change in the screening environment, on the 
timescale of the inverse bandwidth, and an associated reduc¬ 
tion in the size of the Mott gap. The low-energy bosonic 
excitations in the photo-doped Mott insulator also open up 
new relaxation channels, which can substantially reduce the 
doublon-hole recombination or production and hence the ther- 
malization time of the system. 

In strongly correlated metallic systems, the effect of a field 
pulse depends on the pulse frequency and amplitude. In these 
systems, the low-energy quasi-particle band contributes to 
the screening. If the coherent quasi-particles are destroyed 
by a weak low-energy pulse, low-energy screening processes 
are eliminated, and the system can be switched into a more 
strongly interacting, hot insulating state. If the pulse intensity 
is large, or the pulse frequency is of the order of the split¬ 
ting between the Hubbard side-bands, the resulting state is 
essentially a photo-doped narrow-gap insulator in which the 
screening from low-energy excitations within the photo-doped 
Hubbard bands may overcompensate the loss of screening as¬ 
sociated with the vanishing of the quasi-particle peak. Even 
though the effective on-site interaction in such a photo-doped 
narrow gap insulator can be less than in the initial metallic 
equilibrium state, the re-emergence of a quasi-particle peak 


appears to be prevented by heating and strong doublon/hole 
scattering. Within the timescale accessible in our simulations 
(about 30 inverse hoppings) we were only able to observe 
the emergence of a quasi-particle peak in systems coupled to 
a heat bath with broad energy spectrum, or in quenches to 
on-site interactions which are substantially lower than what 
can be achieved by photo-doping in the regime of weak-to- 
intermediate V. 

For a more accurate simulation in the metal-insulator 
crossover regime and the study of screening induced transient 
states near the first order metal-insulator transition, a more 
accurate solver for the electron-boson impurity problem is 
needed. The extension of our NCA + weak coupling based 
solver to one crossing diagrams, or the implementation of the 
real-time OCA + Lang-Firsov scheme would seem the obvi¬ 
ous next steps on the methodological side. 

Our nonequilibrium EDMFT framework is an important 
step towards the development of an ab-initio scheme for out- 
of-equilibrium strongly correlated materials. The screening 
of the Coulomb interaction is a crucially important effect 
in solids, and simple schemes, such as a straight-forward 
nonequilibrium extension of the widely used LDA+DMFT 
framework, will fail to describe changes in the screening en¬ 
vironment induced for example by a laser pulse. Ab-initio 
schemes such as the GW+DMFT methodp^ which are built 
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on top of an EDMFT framework and compute the (partially) 
screened Coulomb interaction in a self-consisten manner, ap¬ 
pear to be a promising route forward. 
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the bosonic propagator W can be expressed a; 


m 


Wi at ')-~2 6HZ) 


= sz * wqM * i 
SU{t\,t 2 ) SU~ x {t',t) z 


(Al) 


where we have used the chain rule and the relations U *U 1 = 
Sc, SU. *U~ X = — U * 5U~ X , which imply 


5U = -U * 5U~ X *U , 

6U~ x (t',t) =- w (*i.* , M*»*2 ). (A 2 ) 


Appendix A: Bosonic propagator from charge-charge 
correlations 

In this appendix we derive the relation between the bosonic 
propagator W and density-density correlator \ defined in 
Eq. {FT} for the nonequilibrium case. Using the action <0 


On the other hand, we may can express the partition func¬ 
tion using Eqs. (JtJi and (Jsj. From this, we obtain - = 

— |Ximp + \Sl~ x and finally arrive at 

W imp =U-U*x imp*W- (A3) 
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